/*****************************************************************************/
/*
/*      Algorithm Name:         ECOH (Elliptic Curve Only Hash)
/*      Principal Submitter:    Daniel R. L. Brown
/*
/*		ecoh.c					This file implements the major cryptographic
/*                              operations for ECOH
/*
/*      Date:                   23 October 2008
/*
/*
/*      COMMAND LINES:
/*
/*	GNU GCC
/*      gcc 64-bit:     gcc -o ecoh.exe -m64 -D "_x64_" -O3 ecoh.param.c ecoh.c ecoh.main.c
/*      gcc 32-bit:     gcc -o ecoh.exe -O3 ecoh.param.c ecoh.c ecoh.main.c
/*		
/*	MICROSOFT VISUAL STUDIO
/*      cl 64-bit		cl /Feecoh /O2 /GL -D "_x64_" ecoh.main.c ecoh.param.c ecoh.c
/*		cl 32-bit		cl /Feecoh /O2 /GL ecoh.main.c ecoh.param.c ecoh.c
/*						
/*****************************************************************************/
#include "ecoh.h"
/*****************************************************************************/
/*
/*	external variables required for this file
/*
/*****************************************************************************/
extern curve	sect283r1;
extern curve	sect409r1;
extern curve	sect571r1;
extern uchar	TBIT[];
extern uchar	BMASK[];
extern uint		BIT[];
extern uint		F2X2[];
extern uint		F2X22[];
extern uint		HT283[][SECT283_SIZE];
extern uint		HT409[][SECT409_SIZE];
extern uint		HT571[][SECT571_SIZE];
extern uint		W_MASK[];
/*****************************************************************************/
/*
/*	Some MACROS for setting and adding counters.  These are provided to make
/*  the code more readable as opposed to the purpose of optimization
/*
/*****************************************************************************/
#ifdef _x64_
#define SET(OP,UI,LEN)		{	(OP)[0]=(UI); if((LEN)==2){	OP[1]=0;	} }
#define SETDL(OP,DL,LEN)	SET((OP),(DL),(LEN))
#define SETCTR(OP,CTR,LEN)	{	(OP)[0]=(CTR)[0]; if((LEN)==2){	(OP)[1]=(CTR)[1];	} }
#define ADD(OP,UI,LEN)		{\
	if((((OP)[0])+=(UI))<(UI)){\
		if((LEN)==2){ if((++((OP)[1]))==0)	return FAIL;	}\
		else return FAIL;\
	}\
}
#define ADDDL(OP,DL,LEN)	( ADD((OP),(DL),(LEN))) 
#define INCCTR(CTR, LEN)	{ \
	if(++((CTR)[0])==0){	\
		if((LEN)==2){\
			if(++((CTR)[1])==0)	return BAD_HASHBITLEN; \
		}\
		else	return BAD_HASHBITLEN;\
	}\
}
#else
#define SET(OP,UI,LEN)		{	(OP)[0]=(UI);	(OP)[1]=0;	if((LEN)==4) { (OP)[2]=0; (OP)[3]=0;	} }
#define SETDL(OP,DL,LEN)	{	(OP)[0]=(uint)((DL)&0xFFFFFFFF);	\
	(OP)[1]=(uint)((DL)>>32);	\
	if((LEN)==4)	{ (OP)[2]=0; (OP)[3]=0;	} \
}
#define SETCTR(OP,CTR,LEN)	{	SETDL((OP), (CTR)[0], 2);	if((LEN)==4)	SETDL(&((OP)[2]), ((CTR)[1]), 2);	}
#define ADD(OP,UI,LEN)		{\
	if((((OP)[0])+=UI)<(UI)){\
		if((++((OP)[1]))==0){\
			if((LEN)==4){\
				if((++((OP)[2]))==0){\
					if((++((OP)[3]))==0)	return FAIL;\
				}}}\
		else	return FAIL;\
	}}\

#define ADDDL(OP,DL,LEN)	{  \
	uint cx = 0;\
	if((((OP)[0])+=(uint)((DL)&0xFFFFFFFF))<(uint)((DL)&0xFFFFFFFF)){\
		cx = 1;\
		if((((OP)[1])+=cx))  cx = 0;\
	}\
	if((((OP)[1])+=(uint)((DL)>>32))<(uint)((DL)>>32))	++cx;\
	if((LEN)==4){\
		if((((OP)[2])+=cx)<cx){\
			if((++((OP)[3]))==0)	return FAIL;\
		}\
	}\
	else \
		if (cx)	return FAIL;\
}
#define INCCTR(CTR, LEN)	{ \
	if(++((CTR)[0])==0){	\
		if((LEN)==4){\
			if(++((CTR)[1])==0)	return BAD_HASHBITLEN; \
		}\
		else	return BAD_HASHBITLEN;\
	}\
}
#endif

#ifdef _DEBUG
void	PrintInt2(uint sz, uint *op)
{	
	uint i;
	for(i=sz-1;i!=-1;i--){	printf(SZ_PRINTINT, op[i]);	}
}
#endif
/*****************************************************************************/
/*
/*	START OF f_{2^m} FUNCTIONS
/*
/*****************************************************************************/
/*
/*	IsInt()				checks to see if int array is = a single int, 
/*
/*  Input Values:   sz      length of the uint array op
/*
/*                  op		pointer to a unit array to evaluate
/*
/*                  ui      uint value to check to see if op value when 
/*                          viewed as a sz-limb value is equal to ui
/*
/*  Return Values	0		on false 
/*
/*                  1		on true
/*
/*****************************************************************************/
int		IsInt(uint sz, uint *op, uint ui)
{
	uint    i;
	if(op[0]!=ui)	return 0;
	for(i=sz-1;i>0;i--){	if(op[i])	return 0;	}
	return 1;
}/****************************************************************************/
/*
/*	GetBitLength()		returns the bit length floor(log_2(op))
/*
/*  Input Values:	sz		length of the uint array op
/*
/*                  op		pointer to a uint array evaluate
/*
/* Return Values   floor(log_2(op))
/*
/*****************************************************************************/
uint	GetBitLength(uint sz, uint *op)
{
	sint	i, d;
	uint	n;
	if(sz == 0)		return 0;
	d = sz-1;
	while(op[d]==0){	--d;	}
	if(d<0)			return 0;	
	n = op[d];
	i = 0;
	while(n){	n>>=1;	++i;	}
	return (d*(LIMB_BIT_SIZE)) + i;
}/****************************************************************************/
/*
/*	PolyShiftRight()	This function shifts a polynomial to the right 1 bit
/*						in F_2[z]
/*
/*	Input Values:	rop		pointer to an uint array to hold the result
/*
/*					sz		length of the uint array op
/*
/*					op		pointer to a uint array to shift to the right
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void	PolyShiftRight(uint *rop, uint sz, uint *op)
{
	uint	t1, t2;
	--sz;
	t2 = 0;
	do{
		if(op[sz]&0x01)	t1 = TOPBIT;
		else			t1 = 0;
		rop[sz] = (op[sz]>>1)|t2;
		t2 = t1;
	}while(--sz!=(-1));

}/****************************************************************************/
/*
/*	PolyAdd()		This function add's to polynomials of the same length in 
/*					F_2[z], rop = op1 + op2 in F_2[z]
/*
/*	Input Values:	rop		pointer to an uint array for the result
/*
/*					sz		length of the arrays op1, op2
/*
/*					op1		pointer to an uint arrary of first operand
/*
/*					op2		pointer to an uint array of second operand
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void	PolyAdd(uint *rop, uint sz, uint *op1, uint *op2)
{
	uint i;
	i = -sz;	op1 -=i;	op2 -=i;	rop -=i;
	do{
		rop[i] = op1[i]^op2[i];
	}while(++i!=0);
}/****************************************************************************/
/*
/*	PolyMulZ()		This function performs a polynomial multiplication 
/*					of op(z) in F_{2}[z] by z^deg.  The result is 
/*
/*					rop(z) =  z^{deg}*op(z).
/*
/*	RESTRICTIONS:	rop must be properly allocated.
/*
/*	Input Values:	rop		pointer to an uint array to hold the result
/*
/*					sz		length of the op array
/*
/*					op		pointer to an uint array to multiply
/*
/*					deg		the degree of z by which to multiply op,
/*							(or shift to the left amount) 
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void	PolyMulZ(uint *rop, uint sz, uint *op, uint deg)
{
	uint			big, r, l, size, n;
	uint			t;

	size = sz + (deg/LIMB_BIT_SIZE);		big = size - sz;
	l = deg - big*LIMB_BIT_SIZE;
	if(l==0){	
		memcpy(rop+big, op, sizeof(uint)*sz);	
		--big;
		while(big!=-1){	rop[big--] = 0;	}
		return;
	}
	r = LIMB_BIT_SIZE - l;
	--sz;
	n = 0;
	rop[sz+big+1] = 0;
	while(sz!=-1){
		t = op[sz];		rop[sz+big] = (t<<l);
		t>>=r;			rop[sz+big+1]|= t;
		--sz;
	}
	--big;
	while(big!=-1){	rop[big--] = 0;	}
}/****************************************************************************/
/*
/*	PolyModT283()	Performs modular reduction by the polynomial reduce
/*					modulo f(z) = z^283 + z^12 + z^7 + z^5 + 1   
/*	
/*					See [1], Algorithm 2.43, page 56
/*
/*	RESTRICTIONS:	deg(op(z)) < 565
/*
/*	Input Values:	rop		pointer to an uint array to hold the result
/*
/*					op		point to an uint array to reduce
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void	PolyModT283(uint *rop, uint *op)
{
	uint j, h1;
#ifdef _x64_
	for(j=8;j>4;j--){
		h1 = op[j];
		op[j-5]^=(h1<<37)^(h1<<42)^(h1<<44)^(h1<<49);
		op[j-4]^=(h1>>27)^(h1>>22)^(h1>>20)^(h1>>15);
	}
	h1 = op[4]>>27;
	op[0] ^= h1^(h1<<5)^(h1<<7)^(h1<<12);
	op[4] &=0x0000000007FFFFFF;
#else
	for(j=17;j>8;j--){
		h1 = op[j];
		op[j-9]^=(h1<<5)^(h1<<10)^(h1<<12)^(h1<<17);
		op[j-8]^=(h1>>27)^(h1>>22)^(h1>>20)^(h1>>15);
	}
	h1 = op[8]>>27;
	op[0] ^= h1^(h1<<5)^(h1<<7)^(h1<<12);
	op[8] &=0x07FFFFFF;
#endif
	memcpy(rop, op, sizeof(uint)*SECT283_SIZE);
}/****************************************************************************/
/*
/*	PolyModT409()	Performs modular reduction by the polynomial reduce
/*					modulo f(z) = z^409 + z^87 + 1   
/*
/*					See [1], Algorithm 2.44, page 56
/*
/*	RESTRICTIONS:	deg(op(z))<817
/*
/*	Input Values:	rop		pointer to an uint array to hold the result
/*
/*					op		point to an uint array to reduce
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void	PolyModT409(uint *rop, uint *op)
{
	uint j, h1;
#ifdef _x64_
	for(j=12;j>6;j--){
		h1 = op[j];						op[j-7]^= (h1<<39);		
		op[j-6]^=(h1<<62)^(h1>>25);		op[j-5]^=(h1>>2);
	}
	h1 = op[6]>>25;						op[0] ^= h1;
	op[1] ^=(h1<<23);					op[6] &=0x0000000001FFFFFF;
#else
	for(j=25;j>12;j--){
		h1 = op[j];						op[j-13]^=(h1<<7);
		op[j-12]^=(h1>>25);				op[j-11]^=(h1<<30);
		op[j-10]^=(h1>>2);
	}
	h1 = op[12]>>25;					op[0] ^= h1;
	op[2] ^= h1<<23;					op[12] &=0x01FFFFFF;
#endif
	memcpy(rop, op, sizeof(uint)*SECT409_SIZE);
}/****************************************************************************/
/*
/*	PolyModT571()	Performs modular reduction by the polynomial reduce
/*					modulo f(x) = x^571 + x^10 + x^5 + x^2 + 1
/*	
/*					See [1], Algorithm 2.45, page 56
/*
/*	RESTRICTIONS:	deg(op(z))<1141
/*
/*	Input Values:	rop		pointer to an uint array to hold the result
/*
/*					op		point to an uint array to reduce
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void	PolyModT571(uint *rop, uint *op)
{
	uint j, h1;
#ifdef _x64_
	for(j=17;j>8;j--){
		h1 = op[j];
		op[j-9]^=(h1<<5)^(h1<<7)^(h1<<10)^(h1<<15);
		op[j-8]^=(h1>>59)^(h1>>57)^(h1>>54)^(h1>>49);
	}
	h1 = op[8]>>59;
	op[0] ^= h1^(h1<<2)^(h1<<5)^(h1<<10);
	op[8] &=0x07FFFFFFFFFFFFFF;
#else
	for(j=35;j>17;j--){
		h1 = op[j];
		op[j-18]^=(h1<<5)^(h1<<7)^(h1<<10)^(h1<<15);
		op[j-17]^=(h1>>27)^(h1>>25)^(h1>>22)^(h1>>17);
	}
	h1 = op[17]>>27;
	op[0] ^= h1^(h1<<2)^(h1<<5)^(h1<<10);
	op[17] &=0x07FFFFFF;
#endif
	memcpy(rop, op, sizeof(uint)*SECT571_SIZE);
}/****************************************************************************/
/*
/*	PolyMulMod()	modular polynomial multiplication in F_2[z]/(f(z))
/*					rop(z) = op1(z)*op2(z) in F_2[z]/(f(z)), where z is implied
/*					by sz the length of operands 1 and 2
/*
/*					performs right-to-left comb method, see [1], Algorithm
/*					2.34, p 49
/*
/*	RESTRICTIONS:	op1(z), op2(z) < f(z)
/*
/*	Input Values:	rop		pointer to an uint array to hold the result
/*
/*					sz		length of the operands op1 and op2
/*		
/*					op1		pointer to a uint array of first operand
/*
/*					op2		pointer to a uint array of first operand
/*
/*	Return Values:	void	
/*
/*****************************************************************************/
void	PolyMulMod(uint *rop, uint sz, uint *op1, uint *op2)
{
	uint	k, j, d[3*MAX_SIZE + 2], *out;
	memset(d, 0, sizeof(uint)*(3*MAX_SIZE+ 2));
	out = d + MAX_SIZE + 2;
	memcpy(d, op2, sizeof(uint)*sz);
	for(k=0;k<LIMB_BIT_SIZE;k++){
		for(j=0;j<sz;j++){
			if(BIT[k]&op1[j])	{	PolyAdd(&out[j], sz+1, &out[j], d);	}
		}
		if(k!=(LIMB_BIT_SIZE-1)){	PolyMulZ(d, sz+1, d, 1);	}
	}
	switch(sz)
	{
	case SECT283_SIZE:		PolyModT283(rop, out);		break;
	case SECT409_SIZE:		PolyModT409(rop, out);		break;
	case SECT571_SIZE:		PolyModT571(rop, out);		break;
	default:											break;
	}
}/****************************************************************************/
/*
/*	PolySqrMod()	modular polynomial square in F_2[z]/(f(z)), where
/*					f(z) is implied by sz the length of op(z),
/*
/*					rop(z) = op(z)^2 (mod f(z))
/*
/*					performs a modified version of [1], Algorithm 2.39, 
/*					p. 53
/*
/*	RESTRICTIONS:	op(z) < f(z)
/*
/*	Input Values:	rop		pointer to an uint array to hold the result
/*
/*					sz		length of the operand op
/*
/*					op		pointer to an uint array to square
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void	PolySqrMod(uint *rop, uint sz, uint *op)
{
	uint	j, out[2*MAX_SIZE];
	/*************************************************************************/
	/*
	/*	square polynomial using the F2X2 look-up table
	/*
	/*************************************************************************/
	for(j=0;j<sz;j++){
#ifdef _x64_
		out[2*j] = F2X2[op[j]&0xFF]|(F2X22[(op[j]&0xFF00)>>8]);;
		out[2*j] |= F2X22[(op[j]&0xFF0000)>>16]<<16|(F2X22[(op[j]&0xFF000000)>>24]<<32);
		out[2*j+1] = F2X2[(op[j]&0xFF00000000)>>32]|(F2X22[(op[j]&0xFF0000000000)>>40]);;
		out[2*j+1] |= F2X22[(op[j]&0xFF000000000000)>>48]<<16|(F2X22[(op[j]&0xFF00000000000000)>>56]<<32);
#else
		out[2*j] = F2X2[op[j]&0xFF]|(F2X22[(op[j]&0xFF00)>>8]);;
		out[2*j+1] = F2X2[(op[j]&0xFF0000)>>16]|(F2X22[(op[j]&0xFF000000)>>24]);
#endif
	}
	switch(sz)
	{
	case SECT283_SIZE:		PolyModT283(rop, out);		break;
	case SECT409_SIZE:		PolyModT409(rop, out);		break;
	case SECT571_SIZE:		PolyModT571(rop, out);		break;
	default:											break;
	}
}/****************************************************************************/
/*
/*	PolyInv()		perform a field inversion (expensive) in F_2[z]/(f(z))
/* 
/*					Performs inversion in F_{2^m} using extended Euclidean,
/*					[1], Algorithm 2.48, p.58
/*
/*	RESTRICTIONS:	op(z) < f(z)
/*	
/*	Input Values:	rop		pointer to an uint array to hold the result
/*
/*					sz		the length of the input operands op and fz
/*
/*					op		pointer to an uint array to invert mod fz
/*
/*					fz		pointer to an uint array to interpret as the
/*							irreducible binary polynomial
/*	
/*	Return Values:	void
/*
/*****************************************************************************/
void	PolyInv(uint *rop, uint sz, uint *op, uint *fx)
{
	uint	scratch[8*(MAX_SIZE+1)], *u, *v, *g1, *g2, *t;
	sint		j;
	t = scratch;
	if(op[0]==0)	{
		if(IsInt(sz, op, 0)){
			memset(rop, 0, sizeof(uint)*sz);		return;
		}
	}
	memset(t, 0, sizeof(uint)*8*(MAX_SIZE+1));
	u = t + 2*(MAX_SIZE+1);				v = u + 2*(MAX_SIZE+1);
	g1 = v + (MAX_SIZE+1);				g2 = g1 + (MAX_SIZE+1);
	memcpy(u, op, sizeof(uint)*sz);		memcpy(v, fx, sizeof(uint)*sz);
	g1[0] = 1;
	while(!IsInt(sz, u, 1)){
		j = GetBitLength(sz, u) - GetBitLength(sz, v);
		if(j<0){
			t = u;		u = v;		v = t;
			t = g1;		g1 = g2;	g2 = t;
			j = -j;		t = scratch;
		}
		PolyMulZ(t, sz, v, j);		PolyAdd(u, sz, u, t);
		PolyMulZ(t, sz, g2, j);		PolyAdd(g1, sz, g1, t);
	};	
	memcpy(rop, g1, sizeof(uint)*sz);
}
/*	END OF f_{2^m} FUNCTIONS
/*
/*****************************************************************************/
/*****************************************************************************/
/*
/*	START OF ECC FUNCTIONS
/*
/*****************************************************************************/
/*
/*	PointDouble()	This function doubles a point represented in
/*					polynomial basis on a curve See [1] Group Laws for 
/*					non-super singular curves p. 81, Point Doubling
/*
/*					rop = op + op on T
/*
/*	Input Values:	rop		pointer to a point to hold result
/*
/*					op		pointer to a point to double
/*
/*					T		pointer to a curve on which to operate
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void PointDouble(point *rop, point *op, curve *T)
{
	uint	l[2*MAX_SIZE], l2[2*MAX_SIZE], px[MAX_SIZE];
	
	if(op->x[0]==0){
		if(IsInt(T->sz, op->x, 0)){
			memset(rop->x, 0, sizeof(uint)*T->sz);
			memset(rop->y, 0, sizeof(uint)*T->sz);
			return;
		}
	}
	PolyInv(l, T->sz, op->x, T->f);		/*	l = x^-1					*/
	PolyMulMod(l, T->sz, l, op->y);		/*	l = y/x in f2m				*/
	PolyAdd(l, T->sz, l, op->x);		/*	l = x + y/x in f2m			*/
	PolySqrMod(l2, T->sz, l);			/*	l2 = l^2 in f2m				*/
	l[0]^=T->a;							/*	l = l + 1 in f2m			*/
	PolyAdd(px, T->sz, l2, l);			/*  px = l^2 + l + a in f2m     */
	PolySqrMod(l2, T->sz, op->x);		/*	l2 = x^2					*/
	PolyMulMod(l, T->sz, l, px);		/*	l = px(l + 1) 				*/
	PolyAdd(rop->y, T->sz, l2, l);		/* py = x^2 + px(l + 1) in f2m  */
	memcpy(rop->x, px, sizeof(uint)*T->sz);
}/****************************************************************************/
/*
/*	PointAdd()		add one point to another on the curve 
/*					See [1] Group Laws for non-super singular curves
/*					p. 81, Point Addition
/*
/*					rop = op1 + op2 on T
/*
/*	Input Values:	rop		pointer to a point to hold the result
/*
/*					op1		pointer to a first point over which to sum
/*
/*					op2		pointer to a second point over which to sum
/*
/*					T		pointer to a curve on which to operate
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void PointAdd(point *rop, point *op1, point *op2, curve *T)
{
	uint	l[2*MAX_SIZE], l2[2*MAX_SIZE], px[MAX_SIZE], py[MAX_SIZE];
	uint	opx0, opy0;

	opx0 = IsInt(T->sz, op1->x, 0);	opy0 = IsInt(T->sz, op1->y, 0);
	if(opx0 && opy0)	{	/*	rop = O + op2	*/
		memcpy(rop->x, op2->x, sizeof(uint)*T->sz);
		memcpy(rop->y, op2->y, sizeof(uint)*T->sz);
		return;
	}
	opx0 = IsInt(T->sz, op2->x, 0);	opy0 = IsInt(T->sz, op2->y, 0);
	if(opx0 && opy0)	{	/*	rop = op1 + O	*/
		memcpy(rop->x, op1->x, sizeof(uint)*T->sz);
		memcpy(rop->y, op1->y, sizeof(uint)*T->sz);
		return;
	}/*	rop = op1 + op1	or rop = op1 + (-op1) = O	*/
	if(memcmp(op1->x, op2->x, sizeof(uint)*T->sz)==0){
		if(memcmp(op1->y, op2->y, sizeof(uint)*T->sz)==0){
			PointDouble(rop, op1, T);		return;
		}
		else{
			memset(rop->x, 0, sizeof(uint)*T->sz);
			memset(rop->y, 0, sizeof(uint)*T->sz);
			return;
		}
	}/*	otherwise rop = op1 + op2	*/
	PolyAdd(l2, T->sz, op1->x, op2->x);	/*	l2 = x1 + x2					*/
	PolyInv(l, T->sz, l2, T->f);		/*	l = 1/(x1 + x2)					*/
	PolyAdd(l2, T->sz, op1->y, op2->y);	/*	l2 = y1 + y2					*/
	PolyMulMod(l, T->sz, l, l2);		/*	l = (y1 + y2)/(x1 + x2)			*/
	PolySqrMod(l2, T->sz, l);			/*	l2 = l^2						*/
	PolyAdd(px, T->sz, l2, l);			/*	px = l^2 + l					*/
	PolyAdd(px, T->sz, px, op1->x);		/*	px = l^2 + l + x1				*/
	PolyAdd(px, T->sz, px, op2->x);		/*	px = l^2 + l + x1 + x2			*/
	px[0]^=T->a;						/*	px = l^2 + l + x1 + x2 + a		*/
	PolyAdd(l2, T->sz, px, op1->x);		/*	l2 = x3 + x1					*/
	PolyMulMod(l2, T->sz, l2, l);		/*	l2 = l(x3 + x1)					*/
	PolyAdd(py, T->sz, l2, px);			/*	py = l(x3 + x1) + x3			*/
	PolyAdd(py, T->sz, py, op1->y);		/*	py = l(x3 + x1) + x3 + y1		*/
	memcpy(rop->x, px, sizeof(uint)*T->sz);	
	memcpy(rop->y, py, sizeof(uint)*T->sz);	
}/****************************************************************************/
/*
/*	BasePointMultiply()		This function performs a base point scalar
/*							multiply. See [1], Algorithm 3.41, p. 104.
/*
/*							rop = k*T->G
/*
/*	Input Values:	rop		pointer to point in which to hold the result
/*
/*					sz		length of the integer array k
/*
/*					k		pointer to an uint array integer
/*			
/*					T		pointer to a curve over which to operate
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void	BasePointMultiply(point *rop, uint sz, uint *k, curve *T)
{
	uint	i, j, n;
	point	B, A;

	n = GetBitLength(sz, k);					n = (n + 3) >> 2;
	memset(A.x, 0, MAX_SIZE*sizeof(uint));	memset(A.y, 0, MAX_SIZE*sizeof(uint));
	memset(B.x, 0, MAX_SIZE*sizeof(uint));	memset(B.y, 0, MAX_SIZE*sizeof(uint));

	for(j=15;j>0;j--){
		for(i=0;i<n;i++){
			if(((k[i/LIMB_W_LEN]&W_MASK[i%LIMB_W_LEN])>>(4*(i%LIMB_W_LEN)))==j){
				PointAdd(&B, &B, &T->G[i], T);		/*	B = B + 2^{4*i}G	*/
			}
		}
		PointAdd(&A, &A, &B, T);							
	}
	memcpy(rop->x, A.x, sizeof(uint)*T->sz);	
	memcpy(rop->y, A.y, sizeof(uint)*T->sz);	
}/****************************************************************************/
/*
/*	SolveQuadratic()			Solves the quadratic equation for rop
/*								rop^2 + rop = op
/*
/*								Uses the basic version algorithm in [1], 
/*								Algorithm 3.85, p.133
/*
/*	RESTRICTIONS:				trace(op) = 0, guarenteeing a solution
/*								op gets destroyed
/*
/*	Input Values:	rop		pointer to an array of uint to hold result
/*
/*					sz		length of op array
/*
/*					op		pointer to an array of uint for which to solve
/*							the quadratic
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void SolveQuadratic(uint *rop, uint sz, uint *op)
{
	uint			i, g[MAX_SIZE], b, qsz, *H;
	memset(g, 0, sizeof(uint)*sz);
	if(sz==SECT283_SIZE)		{	qsz = 141;	H = (uint*)HT283;	}
	else if(sz==SECT409_SIZE)	{	qsz = 204;	H = (uint*)HT409;	}
	else						{	qsz = 285;	H = (uint*)HT571;	}
	for(i=qsz;i>0;i--){
		if(GETBIT(op, 2*i)){
			b = BIT[i%LIMB_BIT_SIZE];
			op[i/LIMB_BIT_SIZE]^= b;
			g[i/LIMB_BIT_SIZE]^= b;
		}
	}
	for(i=0;i<qsz;i++){
		if(GETBIT(op, 2*i + 1)){	
			PolyAdd(g, sz, g, &H[i*sz]);
		}
	}
	memcpy(rop, g, sizeof(uint)*sz);

}/*****************************************************************************/
/*
/*	PointDecompress()	This function attempts to decompress a point. It uses 
/*						a hard coded bit to determine which point on the curve 
/*						it is mapped to.
/*
/*						See [2], Octet-String-to-Elliptic-Curve-Point Conversion
/*						\S 2.3.4, p.11.
/*
/*	RESTRICTIONS:		trace(rop->x) = 1
/*
/*	Input Values:		rop		pointer to input/output of a point to decompress
/*
/*						T		pointer to a curve over which to operate.
/*
/*	Return Values:		0		On Success
/*
/*						-1		On Failure (as UINT)
/*
/*****************************************************************************/
uint	PointDecompress(point *rop, curve *T)
{
	uint	yp, a[2*MAX_SIZE], b[MAX_SIZE], loop = 1;
	/*	compute a = x + T->a + T->b*x^{-2} in F_{2^m}						*/
	PolyInv(b, T->sz, rop->x, T->f);			/*	a = x^-1				*/
	PolySqrMod(a, T->sz, b);					/*	b = x^-2				*/
	PolyMulMod(a, T->sz, a, T->b);				/*	a = T->b*(x^-2)			*/
	switch(T->sz)	{
	case SECT283_SIZE:		if(TRACET283(a))	return -1;	
		yp = GETBIT(rop->x, YP_BIT_T283);		break;
	case SECT409_SIZE:		if(TRACET409(a))	return -1;
		yp = GETBIT(rop->x, YP_BIT_T409);		break;
	case SECT571_SIZE:		if(TRACET571(a))	return -1;
		yp = GETBIT(rop->x, YP_BIT_T571);		break;
	default:
		return -1;
	}
	PolyAdd(a, T->sz, a, rop->x);			/*	a = x + T->b*(x^-2)			*/
	a[0]^=T->a;								/*	a = x + T->a + T->b*(x^-2)	*/
	SolveQuadratic(b, T->sz, a);			/*	solve quadratic z^2 + z = a	*/
	if(yp!=(b[0]&0x01))	{		b[0]^=0x01;	}
	PolyMulMod(a, T->sz, b, rop->x);
	memcpy(rop->y, a, sizeof(uint)*T->sz);
	return 0;
}
/*
/*	END OF ECC FUNCTIONS
/*
/*****************************************************************************/
/*****************************************************************************/
/*
/*	START OF ECOH FUNCTIONS AND AUXILARY FUNCTIONS
/*
/*****************************************************************************/
/*
/*	Unload()			This function unloads a uint array as a BitSequence
/*						according to the ECOH specification [3], although 
/*						perhaps a bit clumsuly.
/*
/*	Input Values:	rop		pointer to a uchar array to hold result
/*
/*					sz		length of the uchar result rop to return
/*
/*					op		pointer to an array of uint from which to copy
/*
/*	Return Values:	void
/*
/*****************************************************************************/
void Unload(uchar *rop, uint sz, uint *op)
{
    uint	i;
	uchar	bval;
	memcpy(rop, (uchar*)op, sz);
	for(i=0;i<(sz/2);i++){
		bval = rop[i];	rop[i] = rop[sz-(i+1)];	rop[sz-(i+1)] = bval;
	}
}
/*****************************************************************************/
/*
/*	FormatData()		This function limb reverses a multi-limb element for
/*						the purpose of formatting prior to the internal hash
/*						transformation, according to ECOH Specification [3]
/*
/*	Input Values:	rop		pointer to an array of uint to format
/*
/*					sz		length of the array rop to format
/*
/*****************************************************************************/
void FormatData(uint *rop, uint sz)
{
	uint	value, ctr, i = 0;

	ctr = sz;
	
	while((--ctr)>(i)){
		value = rop[i];	rop[i++]=rop[ctr];	rop[ctr] = value;
	}
	while( sz-- ){
        value = *rop;
#ifdef _x64_
        value = ((value&0xFF00FF00FF00FF00L)>>8)|((value&0x00FF00FF00FF00FFL)<<8);
        value =((value&0x0000FFFF0000FFFF)<<16)|((value&0xFFFF0000FFFF0000)>>16);
		*rop++= (value>>32)|(value<<32);
#else
		value = ((value&0xFF00FF00L)>>8)|((value&0x00FF00FFL)<<8);
        *rop++ =(value<<16)|(value>>16);
#endif
     }
}/*****************************************************************************/
/*
/*	HashReturn Init(hashState *state, int hashbitlen);
/*
/*		Initializes a hashState with the intended hash length of this particular 
/*		instantiation.  Additionally, any data independent setup is performed, 
/*		conforming to the SHA3 API [4].
/*
/*	Parameters:
/*		
/*		state:		a structure that holds the hashState information
/*	
/*		hashbitlen: an integer value that indicates the length of the hash 
/*					output in bits.
/*
/*	Returns:
/*		
/*		SUCCESS				on success
/*	
/*		BAD_HASHBITLEN		hashbitlen is invalid (e.g., unknown value)
/*
/*****************************************************************************/
HashReturn Init(hashState *state, int hashbitlen)
{
	state->blen = 0;

	state->counter[0] = 0;	state->counter[1] = 0;
	state->mlen[0] = 0;		state->mlen[1] = 0;
	switch(hashbitlen)
	{
	case ECOH224:
		state->type = ECOH224;		state->DATASIZE = ECOH224_DATASIZE;
		state->clen = ECOH224_CLEN;	state->DIGESTSIZE = ECOH224_DIGESTSIZE;	
		state->T = &sect283r1;
		break;
	case ECOH256:
		state->type = ECOH256;		state->DATASIZE = ECOH256_DATASIZE;
		state->clen = ECOH256_CLEN;	state->DIGESTSIZE = ECOH256_DIGESTSIZE;
		state->T = &sect283r1;
		break;
	case ECOH384:
		state->type = ECOH384;		state->DATASIZE = ECOH384_DATASIZE;
		state->clen = ECOH384_CLEN;	state->DIGESTSIZE = ECOH384_DIGESTSIZE;
		state->T = &sect409r1;
		break;
	case ECOH512:
		state->type = ECOH512;		state->DATASIZE = ECOH512_DATASIZE;
		state->clen = ECOH512_CLEN;	state->DIGESTSIZE = ECOH512_DIGESTSIZE;
		state->T = &sect571r1;
		break;
	default:
		return BAD_HASHBITLEN;
	}
	memset(state->N, 0, MAX_DATASIZE*LIMB_SIZE);
	memset(state->Q.x, 0, MAX_SIZE*LIMB_SIZE);
	memset(state->Q.y, 0, MAX_SIZE*LIMB_SIZE);
	memset(state->P.x, 0, MAX_SIZE*LIMB_SIZE);
	memset(state->P.y, 0, MAX_SIZE*LIMB_SIZE);
	return SUCCESS;
}
/*****************************************************************************/
/*
/*	_Transform()		internal hashing of the data according to the ECOH
/*						specification [3].
/*
/*	Input Values:	state	pointer to the hash state on which to operate
/*
/*	Return Values:	SUCCESS				on Success
/*	
/*					BAD_HASHBITLEN		if maximum hash length breached
/*
/*					FAIL				otherwise
/*
/*****************************************************************************/
uint	_Transform(hashState *state)
{
	SET(state->P.x, 1, state->clen);
	SETCTR(&state->P.x[state->clen], state->counter, state->clen);
	while(PointDecompress(&state->P, state->T))	{
		ADD(state->P.x, 2, state->clen);
	}/*		Q = Q + P, up counter, zeroize blen	*/
	PointAdd(&state->Q, &state->Q, &state->P, state->T);	
#ifdef _DEBUG
	{
		printf("P_%d\t(", (unsigned int)state->counter[0]);	
		PrintInt2(state->T->sz, state->P.x);
		printf(",\n   \t ");
		PrintInt2(state->T->sz, state->P.y);
		printf(")\n\nQ_%d\t(", (unsigned int)state->counter[0]);	
		PrintInt2(state->T->sz, state->Q.x);
		printf(",\n   \t ");
		PrintInt2(state->T->sz, state->Q.y);
		printf(")\n\nN\t ");
		PrintInt2(state->DATASIZE/sizeof(uint), state->N);
		printf("\n\n");
	}
#endif
	INCCTR(state->counter, state->clen);
	state->blen = 0;
	return SUCCESS;
}/*****************************************************************************/
/*
/*	Update()			This function is the update routine for ECOH 
/*						Specification [3], according to the SHA3 API [4].
/*
/*	HashReturn Update(	hashState *state, const BitSequence *data, DataLength 
/*						databitlen);
/*
/*	Process the supplied data.
/*
/*	Parameters:
/*
/*		state:		a structure that holds the hashState information
/*
/*		data:		the data to be hashed
/*
/*		databitlen: the length, in bits, of the data to be hashed
/*
/*	Returns:
/*
/*		SUCCESS			on success
/*	
/*		FAIL			on failure
/*	
/*		BAD_HASHBITLEN	on bad hash bit length
/*
/*****************************************************************************/
HashReturn Update(hashState *state, const BitSequence *data, DataLength databitlen)
{
	DataLength	free, used;
	uint		datastart;
	uchar		*ptr;

	if(databitlen&0x7)			{	return BAD_HASHBITLEN;		}
	if(state->type == ECOH512)	{	datastart = DATASTART*2;	}
	else						{	datastart = DATASTART;		}
	/*	add in the bit count	*/
	state->mlen[0]+=(databitlen);
	if(state->mlen[0]<databitlen)	if((++state->mlen[1])==0)	return BAD_HASHBITLEN;	
	databitlen>>=3;
	used = state->blen;
	ptr = ((unsigned char*)(&state->P.x[datastart])) + used;
	if(used){	/*	handle the case where partial data is in the buffer	*/
		free = state->DATASIZE - used;
		if(databitlen>=free){
			memcpy(ptr, data, (uint)free);		/*	copy into the buffer	*/
			databitlen-=free;			databitlen+=free;
			/*	format the data		*/
			FormatData(&state->P.x[datastart], (uint)state->DATASIZE/sizeof(uint));
			/*	compute N value		*/
			PolyAdd(state->N, state->DATASIZE/sizeof(uint), state->N, &state->P.x[datastart]);
			/*	transform			*/
			if(_Transform(state))	return FAIL;
		}
		else{/*	copy remaining data in buffer and exit		*/
			memcpy(ptr, data, (uint)databitlen);	
			state->blen += (uint)databitlen;
			return SUCCESS;
		}
	}/*	loop over the rest of the data, transforming every ECOH_DATASIZE	*/
	while(databitlen>=state->DATASIZE){
		memcpy(&state->P.x[datastart], data, (uint)state->DATASIZE);
		/*	format the data			*/
		FormatData(&state->P.x[datastart], (uint)state->DATASIZE/sizeof(uint));
		/*	compute N value			*/
		PolyAdd(state->N, state->DATASIZE/sizeof(uint), state->N, &state->P.x[datastart]);
		/*	transform				*/
		if(_Transform(state))	return FAIL;
		databitlen-=state->DATASIZE;
		data+=state->DATASIZE;
	}/*	Load any remaining data		*/
	if(databitlen>0){/*	copy any residual data		*/
		memcpy(&state->P.x[datastart], data, (uint)databitlen);
		state->blen = (uint)databitlen;
	}
	return SUCCESS;
}/****************************************************************************/
/*
/*	Final()			Executes the final operation of ECOH Specification [3], 
/*					according to the SHA3 API [4].
/*
/*		HashReturn Final(hashState *state, BitSequence *hashval);
/*
/*		Perform any post processing and output filtering required and return 
/*		the final hash value.
/*
/*	Parameters:
/*
/*		state:		a structure that holds the hashState information
/*
/*		hashval:	the storage for the final hash value to be returned
/*
/*	Returns:
/*
/*		SUCCESS		on success
/*
/*		FAIL		otherwise
/*	
/*****************************************************************************/
HashReturn Final(hashState *state, BitSequence *hashval)
{
    uint		cnt, datastart;
	uchar	    *ptr;

	if(state->type == ECOH512)	{	datastart = DATASTART*2;	}
	else						{	datastart = DATASTART;		}
	cnt = state->blen;
	ptr = ((unsigned char*)(&state->P.x[datastart])) + cnt;
    *ptr++ = 0x80;
    cnt = state->DATASIZE - 1 - cnt;
    memset(ptr, 0, cnt );
	/*	format data				*/
    FormatData(&state->P.x[datastart], (uint)state->DATASIZE/sizeof(uint));
	/*	compute N value			*/
	PolyAdd(state->N, state->DATASIZE/sizeof(uint), state->N, &state->P.x[datastart]);
	/*	transform				*/
	if(_Transform(state))	return FAIL;
	memcpy(&state->P.x[datastart], state->N, state->DATASIZE);
	SET(state->P.x, 1, state->clen);
	SETCTR(&state->P.x[state->clen], state->mlen, state->clen);
	while(PointDecompress(&state->P, state->T))	{
		ADD(state->P.x, 2, state->clen);
	}/*	Q = Q + P, up counter, zeroize blen					*/
	PointAdd(&state->Q, &state->Q, &state->P, state->T);
	/*	Q = Q + [Q.x/2]G where[Q.x/2] interpreted as an integer and divided by
	/*	2, and G is the base point											*/
	PolyShiftRight(state->P.x, state->T->sz, state->Q.x);
	BasePointMultiply(&state->P, state->T->sz, state->P.x, state->T);
	PointAdd(&state->Q, &state->Q, &state->P, state->T);
	PolyShiftRight(state->Q.x, state->T->sz, state->Q.x);
	/*	unload and return				*/
	Unload(hashval, state->DIGESTSIZE, state->Q.x);
	return SUCCESS;
}/****************************************************************************/
/*	Hash()				This function implements an all-in-one call hash 
/*						algorithm of ECOH [3] using the SHA3 API[4].
/*
/*		HashReturn Hash(int hashbitlen, const BitSequence *data,  DataLength 
/*							databitlen, BitSequence *hashval);
/*
/*
/*		Hash the supplied data and provide the resulting hash value. Set return 
/*		code as appropriate.
/*
/*	Parameters:
/*
/*		hashbitlen:			the length in bits of the desired hash value
/*
/*		data:				the data to be hashed
/*
/*		databitlen:			the length, in bits, of the data to be hashed
/*
/*		hashval:			the resulting hash value of the provided data
/*
/*	Returns:
/*
/*		SUCCESS				on success
/*
/*		FAIL				arbitrary failure
/*
/*		BAD_HASHBITLEN		unknown hashbitlen requested
/*
/*****************************************************************************/
HashReturn Hash(int hashbitlen, const BitSequence *data,  DataLength databitlen, BitSequence *hashval)
{
	hashState	state;
	DataLength	blen, dblen, cnt;
	uchar		*ptr;
	uint		datastart;

	if(Init(&state, hashbitlen))	{	return BAD_HASHBITLEN;	}
	if(state.type == ECOH512)	{	datastart = DATASTART*2;	}
	else						{	datastart = DATASTART;		}
	/*	set mlen, and the byte length blen, and the databitlength remainder	*/
	state.mlen[0]=(databitlen);
	blen = (databitlen>>3);
	dblen = (databitlen&7);
	/*	loop over the data, transforming every ECOH_DATASIZE				*/
	while(blen>=state.DATASIZE){
		memcpy(&state.P.x[datastart], data, (uint)state.DATASIZE);
		/*	format the data			*/
		FormatData(&state.P.x[datastart], (uint)state.DATASIZE/sizeof(uint));
		/*	compute N value			*/
		PolyAdd(state.N, state.DATASIZE/sizeof(uint), state.N, &state.P.x[datastart]);
		/*	transform				*/
		if(_Transform(&state))	return FAIL;
		blen-=state.DATASIZE;
		data+=state.DATASIZE;
	}/*	Load any remaining data		*/
	if(blen>0){
		memcpy(&state.P.x[datastart], data, (uint)blen);
		state.blen = blen;		data+=blen;
	}
	cnt = state.blen;
	ptr = ((unsigned char*)(&state.P.x[datastart])) + cnt;
    if(dblen)	{		*ptr++ = ((*data)&BMASK[dblen])|TBIT[dblen];	}
	else		{		*ptr++ = 0x80;									}

    cnt = state.DATASIZE - 1 - cnt;
    memset(ptr, 0, cnt );
	/*	format 2nd to last block	*/
	FormatData(&state.P.x[datastart], (uint)state.DATASIZE/sizeof(uint));
	/*	compute N value				*/
	PolyAdd(state.N, state.DATASIZE/sizeof(uint), state.N, &state.P.x[datastart]);
	/*	transform					*/
	if(_Transform(&state))	return FAIL;
	/*	now load N and do final block	*/
    memcpy(&state.P.x[datastart], state.N, state.DATASIZE);
	SET(state.P.x, 1, state.clen);
	SETCTR(&state.P.x[state.clen], state.mlen, state.clen);
	while(PointDecompress(&state.P, state.T))	{
		ADD(state.P.x, 2, state.clen);
	}/*	Q = Q + P, up counter, zeroize blen			*/
	PointAdd(&state.Q, &state.Q, &state.P, state.T);
#ifdef _DEBUG
	{
		printf("P_%d\t(", (unsigned int)state.counter[0]);	
		PrintInt2(state.T->sz, state.P.x);
		printf(",\n   \t ");
		PrintInt2(state.T->sz, state.P.y);
		printf(")\n\nQ_%d\t(", (unsigned int)state.counter[0]);	
		PrintInt2(state.T->sz, state.Q.x);
		printf(",\n   \t ");
		PrintInt2(state.T->sz, state.Q.y);
		printf(")\n\nN  \t ");
		PrintInt2(state.DATASIZE/sizeof(uint), state.N);
}
#endif
	/*	Q = Q + [Q.x/2]G where[Q.x/2] interpreted as an integer and divided by
	/*	2, and G is the base point											*/
	PolyShiftRight(state.P.x, state.T->sz, state.Q.x);
#ifdef _DEBUG
	{
		printf("\n\nQ_%d.x/2\t ", (unsigned int)state.counter[0]);	
		PrintInt2(state.T->sz, state.P.x);	
	}
#endif
	BasePointMultiply(&state.P, state.T->sz, state.P.x, state.T);
	PointAdd(&state.Q, &state.Q, &state.P, state.T);
#ifdef _DEBUG
	{
		printf("\n\n(Q_%d.x/2)*G\t(", (unsigned int)state.counter[0]);	
		PrintInt2(state.T->sz, state.P.x);
		printf(",\n           \t ");
		PrintInt2(state.T->sz, state.P.y);
		printf(")\n\nQ_%d + (Q_%d.x/2)*G\t(", (unsigned int)state.counter[0], (unsigned int)state.counter[0]);	
		PrintInt2(state.T->sz, state.Q.x);
		printf(",\n                \t ");
		PrintInt2(state.T->sz, state.Q.y);
	}
#endif
	PolyShiftRight(state.Q.x, state.T->sz, state.Q.x);
#ifdef _DEBUG
	{
		printf(")\n\n(Q + (Q_%d.x/2)*G).x/2\t ", (unsigned int)state.counter[0]);	
		PrintInt2(state.T->sz, state.Q.x);
		printf("\n\n");
	}
#endif
	/*	unload and return		*/
	Unload(hashval, state.DIGESTSIZE, state.Q.x);
	return SUCCESS;
}/*	END OF ECOH FUNCTIONS
/*
/*****************************************************************************/

/*****************************************************************************/
/*
/*	[1] Guide to Elliptic Curve Cryptography, Hankerson, D., Menezes, A., 
/*		Vanstone, S., Springer-Verlag, New York, 2004.
/*
/*	[2] Standards for Efficient Cryptography Group, SEC 1, Elliptic Curve 
/*		Cryptography, Working Draft, Version 1.9, Brown, D.L.R., 22 August 2008,
/*		http://www.secg.org/download/aid-773/sec1_1point9.pdf 
/*
/*	[3] ECOH: the Elliptic Curve Only Hash, Brown, D.L.R., 23 October 2008.
/*
/*	[4]	ANSI C Cryptographic Profile for SHA-3 Candidate Algorithm Submission,
/*		Revision 5, 11 February 2008, NIST,
/*		http://csrc.nist.gov/groups/ST/hash/documents/SHA3-C-API.pdf
/*
/*****************************************************************************/
